
include("../src/qmcmary.jl")
using ..qmcmary

using Dates
using LinearAlgebra

BLAS.set_num_threads(1)


"""保存hs场的位型"""
function savehs(hscfg, head, fname)
    fil = open(fname, "a")
    write(fil, head)
    write(fil, "\n")
    Nt = size(hscfg)[1]
    for ni in Base.OneTo(Nt)
        write(fil, string(hscfg[ni, :]))
        write(fil, "\n")
    end
    close(fil)
end


function chai(L, Nt)
    #the = π*rand()-0.5π
    the = [0.72]
    #Nx = 2*L^2
    Nx = L
    nx = zeros(Nx)
    ny = sin.(the)
    nz = cos.(the)
    println(ny)
    println(nz)
    #
    #hk = lattice_hexagonal(ComplexF64, L, -1.0+0.0im; λ=0.0+0.0im)
    #hk = lattice_chain(L, -1.0+0.0im)
    hk = lattice_ionic_chain(ComplexF64, L, -1.0+0.0im, 0.5+0.0im)
    #println(hk)
    #return
    engs = eigvals(hk)
    println(engs)
    ##
    Ui = 4*ones(Nx)
    Ng = 3
    bidxs = Vector{Int}(undef, Ng)
    hscfg = Matrix{Int}(undef, Nt, Nx)
    for bi in Base.OneTo(Nt)
        cfg = 2(round.(rand(Nx)) .- 0.5)
        hscfg[bi, :] .= cfg
    end
    #hscfg = [-1; +1; -1; +1; -1;;]
    ###
    sp = default_splitting3(Nt, hk, Ui; Z2=true)
    sslen, bmats, allbmats, ss = initialize_SS(Nt, Ng, sp, hscfg, nx, ny, nz)
    ttbar = zeros(Nx)
    nup_i = zeros(2*Nx)
    tph = 0.0
    #
    rightnow = Dates.now()
    fname = "hscfg"*string(Dates.hour(rightnow))*string(Dates.minute(rightnow))
    write(fname, string(the)*"\n")
    #
    nsp = 5000
    for _ in Base.OneTo(1000)
        hscfg, bmats, allbmats, ss, gf, ph = dqmc_step(Nt, Ng, sp, hscfg, nx, ny, nz, sslen, bmats, allbmats, ss)
    end
    #
    tapemap = pbpn_calc_tapemap(sp, nx, ny, nz)
    for idx in Base.OneTo(nsp)
        hscfg, bmats, allbmats, ss, gf, ph = dqmc_step(Nt, Ng, sp, hscfg, nx, ny, nz, sslen, bmats, allbmats, ss)
        println(ph)
        tph += ph
        nxbar, nybar, nzbar = meas_grad(ss, allbmats, sp, hscfg, nx, ny, nz; tapemap=tapemap)
        thetabar = @. nybar * cos(the) -  nzbar * sin(the)
        #println(thetabar)
        #println(diag(gf))
        ttbar += real(thetabar)
        nup_i += diag(gf)
        savehs(hscfg, string(idx), fname)
    end
    println(the, " ", ttbar/nsp)
    println(tph/nsp)
    println(nup_i/nsp)
    fil = open("info", "a")
    write(fil, string(rightnow)*"\n")
    write(fil, string(the)*" "*string(ttbar/nsp)*"\n")
    write(fil, string(tph/nsp)*"\n")
    write(fil, string(nup_i/nsp)*"\n")
end


#=
λ=1.0时，最优在A=（0.57， 0.82）B=（-0.57，0.82）
theta = 0.609385
=#
function optimal_theta1(L)
    #
    #the = 0.609385
    Nx = 2*L^2
    #
    the = zeros(Nx)
    for idx in Base.OneTo(Nx)
        println(-1, " ", (idx+1))
        the[idx] = ^(-1, idx+1)
    end
    #the = [0.6040877380745839, -0.6154628641129916, 0.6138892027453763, -0.5591857376108517, 0.5986068660664536, -0.5970444964954694, 0.5931830662789608, -0.6272841430938022]
    the = rand(Nx)#the*0.72
    #
    nx = zeros(Nx)
    ny = sin.(the)
    nz = cos.(the)
    println(the)
    println(ny)
    println(nz)
    #
    hk = lattice_hexagonal(ComplexF64, L, -1.0+0.0im; λ=sqrt(3)+0.0im)
    Ui = 10*ones(Nx)
    Nt = 60
    Ng = 4
    hscfg = Matrix{Int}(undef, Nt, Nx)
    for bi in Base.OneTo(Nt)
        cfg = 2(round.(rand(Nx)) .- 0.5)
        hscfg[bi, :] .= cfg
    end
    sp = default_splitting2(Nt, hk, Ui; Z2=false)
    sslen, bmats, allbmats, ss = initialize_SS(Nt, Ng, sp, hscfg, nx, ny, nz)
    #
    for _ in Base.OneTo(1)
        hscfg, bmats, allbmats, ss, gf, ph = dqmc_step(Nt, Ng, sp, hscfg, nx, ny, nz, sslen, bmats, allbmats, ss)
        println(ph)
    end
    #
    rightnow = Dates.now()
    fname = "hscfg"*string(Dates.hour(rightnow))*string(Dates.minute(rightnow))
    write(fname, string(the)*"\n")
    #
    ttbar = zeros(Nx)
    nup_i = zeros(2*Nx)
    tph = 0.0
    nsp = 100
    for idx in Base.OneTo(nsp)
        hscfg, bmats, allbmats, ss, gf, ph = dqmc_step(Nt, Ng, sp, hscfg, nx, ny, nz, sslen, bmats, allbmats, ss)
        tph += ph
        println(ph)
        #nxbar, nybar, nzbar = meas_grad(ss, allbmats, ehk, hscfg, Ui, nx, ny, nz; Z2=false)
        #thetabar = @. nybar * cos(the) -  nzbar * sin(the)
        #println(thetabar)
        #ttbar += real(thetabar)
        #println(diag(gf))
        nup_i += diag(gf)
        savehs(hscfg, string(idx), fname)
    end
    #println(the, " ", ttbar/nsp)
    println(tph/nsp)
    println(nup_i/nsp)
    fil = open("info", "a")
    write(fil, string(rightnow)*"\n")
    #write(fil, string(the)*" "*string(ttbar/nsp)*"\n")
    write(fil, string(tph/nsp)*"\n")
    write(fil, string(nup_i/nsp)*"\n")
end


#=
在λ=√2 的时候，最优的位置可以取
=#


#chai(3, 40)
#optimal_theta1(6)

function optimal_psi_the(L)
    kag = lattice_kagome(ComplexF64, L, -1.0+0.0im)
    println(eigvals(kag))
    println(eigvecs(kag)[:, 1])
end

chai(1, 3)
#optimal_theta1(2)
#optimal_psi_the(3)
